!dk parmov
      subroutine parmove 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use corgan_com_M 
      use cindex_com_M 
      use numpar_com_M 
      use cprplt_com_M, ONLY:  
      use cophys_com_M 
      use blcom_com_M, ONLY: qom, t_wall, iphd2, wate, bxv, byv, bzv, ex, ey, &
         ez, x, y, z, tsix, tsiy, tsiz, iphead, ijkcell, ijktmp2, ijkctmp, &
         nux, nuy, nuz, etax, etay, etaz, px, py, pz, pxi, peta, pzta, up, vp, &
         wp, link, ico, uptilde, vptilde, wptilde, bxpn, bypn, bzpn, expn, &
         eypn, nptotl, ezpn, vrms 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j1 
      integer , dimension(8) :: iv 
      integer , dimension(4) :: lioinput, lioprint 
      integer :: lpx, n, np, is, inew, jnew, knew, ijknew 
      real :: nuxp, nuyp, nuzp, nu 
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2 
      real(double) :: ixi4, ieta4, ixi5, ieta5, dto2, qomdt, omdtsq, denom, ut&
         , vt, wt, udotb 
      logical :: expnd, nomore 
      character :: name*80 
!-----------------------------------------------
!
!
!
!     *******************************************
!
!     an explicit particle mover for particles in a static
!     magnetic field
!     the equations of motion are solved in natural coordinates
!
!     *********************************************
      lpx = nphist + 1 
      dto2 = 0.5*dt 
!
      iphd2(ijkcell(:ncells)) = 0 
!
      nptotl = 0 
!
    1 continue 
      nomore = .TRUE. 
!
!     check for more particles to process
!
      ncellsp = 0 
      do n = 1, ncells 
         if (iphead(ijkcell(n)) /= nphist) cycle  
         ncellsp = ncellsp + 1 
         ijkctmp(ncellsp) = ijkcell(n) 
         nptotl = nptotl + 1 
!
      end do 
!
!     if there are more, process particles as though there were one
!     in every cell in the mesh
!
      if (ncellsp /= 0) then 
!
         nomore = .FALSE. 
!
!     for leap-frog, initialize particle positions at half-time level
!
         if (ncyc == 1) then 
            do n = 1, ncellsp 
!
               ijk = ijkctmp(n) 
               np = iphead(ijk) 
!
               px(np) = px(np) + up(np)*dto2 
               py(np) = py(np) + vp(np)*dto2 
               pz(np) = pz(np) + wp(np)*dto2 
!
            end do 
!
            call parlocat (ncellsp, ijkctmp, iphead, iwid, jwid, kwid, nsampl, &
               vrms, t_wall, ico, ijktmp2, rmaj, dz, itdim, npart, ibar, jbar, &
               kbar, mgeom, cdlt, sdlt, wate, x, y, z, tsix, tsiy, tsiz, etax, &
               etay, etaz, nux, nuy, nuz, px, py, pz, up, vp, wp, pxi, peta, &
               pzta, nu_len, nu_comm) 
!
         endif 
!
         call watev (ncellsp, ijkctmp, iphead, itdim, pxi, peta, pzta, wate) 
!
!
         call trilin (ncellsp, ijkctmp, itdim, iwid, jwid, kwid, wate, ex, ey, &
            ez, expn, eypn, ezpn) 
!
         call trilin (ncellsp, ijkctmp, itdim, iwid, jwid, kwid, wate, bxv, byv&
            , bzv, bxpn, bypn, bzpn) 
!
!
!
!     **************************************************************
!
!     solve the momentum equation
!
!     ************************************************************
!
         do n = 1, ncellsp 
!
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
!
            is = ico(np) 
            qomdt = qom(is)*dto2 
!
!
            omdtsq = qomdt**2*(bxpn(ijk)**2+bypn(ijk)**2+bzpn(ijk)**2) 
            denom = 1./(1. + omdtsq) 
!
!
!     solve the momentum equation
!
            ut = up(np) + qomdt*expn(ijk) 
            vt = vp(np) + qomdt*eypn(ijk) 
            wt = wp(np) + qomdt*ezpn(ijk) 
!
            udotb = ut*bxpn(ijk) + vt*bypn(ijk) + wt*bzpn(ijk) 
!
            uptilde(ijk) = (ut + qomdt*(vt*bzpn(ijk)-wt*bypn(ijk)+qomdt*udotb*&
               bxpn(ijk)))*denom 
            vptilde(ijk) = (vt + qomdt*(wt*bxpn(ijk)-ut*bzpn(ijk)+qomdt*udotb*&
               bypn(ijk)))*denom 
            wptilde(ijk) = (wt + qomdt*(ut*bypn(ijk)-vt*bxpn(ijk)+qomdt*udotb*&
               bzpn(ijk)))*denom 
!
         end do 
!
!
         do n = 1, ncellsp 
!
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
!
            up(np) = 2.*uptilde(ijk) - up(np) 
            vp(np) = 2.*vptilde(ijk) - vp(np) 
            wp(np) = 2.*wptilde(ijk) - wp(np) 
!
            px(np) = px(np) + up(np)*dt 
            py(np) = py(np) + vp(np)*dt 
            pz(np) = pz(np) + wp(np)*dt 
!
         end do 
!
!
!     ******************************************************************
!
!     calculate new natural coordinates
!
!     ******************************************************************
!
         call parlocat (ncellsp, ijkctmp, iphead, iwid, jwid, kwid, nsampl, &
            vrms, t_wall, ico, ijktmp2, rmaj, dz, itdim, npart, ibar, jbar, &
            kbar, mgeom, cdlt, sdlt, wate, x, y, z, tsix, tsiy, tsiz, etax, &
            etay, etaz, nux, nuy, nuz, px, py, pz, up, vp, wp, pxi, peta, pzta&
            , nu_len, nu_comm) 
!
!
!
         do n = 1, ncellsp 
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
            if (np <= 0) cycle  
!
            inew = int(pxi(np)) 
            jnew = int(peta(np)) 
            knew = int(pzta(np)) 
!
            ijknew = (knew - 1)*kwid + (jnew - 1)*jwid + (inew - 1)*iwid + 1 
!
            iphead(ijk) = link(np) 
            link(np) = iphd2(ijknew) 
            iphd2(ijknew) = np 
!
         end do 
!
!     ******************************************************************
!
!     all particles have been processed
!
!     ******************************************************************
 
         if (.not.nomore) go to 1 
!
      endif 
!
      do n = 1, ncells 
         ijk = ijkcell(n) 
         iphead(ijk) = iphd2(ijk) 
         iphd2(ijk) = 0 
      end do 
!
!
!     nh=mod(ncyc,nhst)
!     xphys(nh)=px(nphist)
!     yphys(nh)=py(nphist)
!     zphys(nh)=pz(nphist)
!     upar(nh)=up(nphist)
!     vpar(nh)=vp(nphist)
!     wpar(nh)=wp(nphist)
!     epar(nh)=upar(nh)**2+vpar(nh)**2+wpar(nh)**2
!     bxpar(nh)=bxp
!     bypar(nh)=byp
!     bzpar(nh)=bzp
!     bmag=sqrt(bxp**2+byp**2+bzp**2)+1.e-20
!     udotb=upar(nh)*bxp+vpar(nh)*byp+wpar(nh)*bzp
!     usq=upar(nh)**2+vpar(nh)**2+wpar(nh)**2
!     mupar(nh)=(usq-udotb**2/bmag**2)/bmag
!
!
      return  
!l    ------------------------------------------> return ----->>>
!
      end subroutine parmove 
